function [Lambda, RCs, Stand] = mySSA_decompose(Input, M)
	% [Lambda, RCs, Stand] = mySSA_decompose(Input, M)
    %
	% SSA分解函数
    %   1. 将时间序列按照SSA算法分解成重建因子。
    %   2. 如果时间序列出现缺失，会报错。
	% 输入变量：
	% 	Input： 输入的时间序列，需要是列向量
	%	M：		分解时用的窗口大小
	% 输出变量：
	%	Lambda：	特征值，列向量
	%	RCs：		重建的信号序列
	%	Stand:		标准化的振幅、截距
	% 作者：sid（cnqdcyq@163.com）
	% 时间：2021/7/28
	TimeRange = length(Input);

	% 检查是否序列完整
	Input = myMissingCheck(Input);

	% 检查是否是列向量
	Input = my_columncheck(Input);
	% 标准化
	[Input, Stand] = subStandardization(Input);
	% Toeplitz Matrix
	covX = xcorr(Input,M-1,'unbiased');
	Ctoep = toeplitz(covX(M:end));
	% Calculate eigenvalue Lambda and eigenvector RHO
	[RHO,Lambda] = eig(Ctoep);
	Lambda = diag(Lambda);
	[Lambda,index] = sort(Lambda,'descend');
	RHO = RHO(:,index);
	% Calculate Principle Components PCs
	N = TimeRange;
	Y=zeros(N-M+1,M);
	for m=1:M    
		Y(:,m) = Input((1:N-M+1)+m-1);
	end
	PC = Y * RHO;
	%% Calculate reconstructed components RCs
	RCs = zeros(N,M);
	for m=1:M
	    buf=PC(:,m)*RHO(:,m)'; % invert projection
	    buf=buf(end:-1:1,:);
	       for n=1:N % anti-diagonal averaging
	          RCs(n,m)=mean( diag(buf,-(N-M+1)+n) );  
	       end
	end
end

function [output, Stand] = subStandardization(input)
	% 子函数，标准化
	% 输入变量：
	%	input：	列向量，待处理的时间序列
	% 输出变量：
	%	output：	列向量，标准化后的时间序列
	% 	Stand：		标准化的方差和截距[Mean, Std]	
	Mean = mean(input);
	output = input - Mean;
	Std = std(output,1);
	output = output ./ Std;
	Stand = [Mean,Std];
end

function output = myMissingCheck(input)
	if any(isnan(input))
		error("有缺失值");
    end
    output = input;
end
